#ifndef Matrix_
#define Matrix_

#include <iostream>

#include "AlogException.h"


template<class T>
class Matrix {
  public:
    Matrix( int r = 0, int c = 0 );
    Matrix( const Matrix<T>& m );
    ~Matrix() { delete [] element; }

    int Rows() const { return rows; }
    int Columns() const { return cols; }

    T& operator()( int i, int j ) const;
    Matrix<T>& operator=( const Matrix<T>& m );
    Matrix<T> operator+() const;
    Matrix<T> operator+( const Matrix<T>& m ) const;
    Matrix<T> operator-() const;
    Matrix<T> operator-( const Matrix<T>& m ) const;
    Matrix<T> operator*( const Matrix<T>& m ) const;
    Matrix<T>& operator+=( const T& x );

    void Output( std::ostream& out ) const;
  private:
    int rows, cols;
    T *element;
};

template<class T>
Matrix<T>::Matrix( int r/*= 0 */, int c/*= 0 */ )
{
    if ( r < 0 || c < 0 ) throw BadInitializers();

    if ( (!r || !c) && (r || c) )
        throw BadInitializers();

    rows = r; cols = c;
    element = new T[r * c];
}

template<class T>
Matrix<T>::Matrix( const Matrix<T>& m )
{
    rows = m.rows; cols = m.cols;
    element = new T[rows * cols];

    for ( int i = 0; i < rows * cols; i++ )
        element[i] = m.element[i];
}

template<class T>
Matrix<T>& Matrix<T>::operator=(const Matrix<T>& m)
{
    if ( this != &m ) {
        delete [] element;
        rows = m.rows; cols = m.cols;
        element = new T[rows * cols];

        for ( int i = 0; i < rows * cols; i++ )
            element[i] = m.element[i];
    }

    return *this;
}

template<class T>
T& Matrix<T>::operator()( int i, int j ) const
{
    if ( i < 1 || i > rows || j < 1
         || j > cols ) throw OutOfBounds();

    return element[(i - 1) * cols + j - 1];
}

template<class T>
Matrix<T> Matrix<T>::operator+( const Matrix<T>& m ) const
{
    if ( rows != m.rows || cols != m.cols )
        throw SizeMismatch();

    Matrix<T> w( rows, cols );

    for ( int i = 0; i < rows * cols; i++ )
        w.element[i] = element[i] + m.element[i];

    return w;
}

template<class T>
Matrix<T> Matrix<T>::operator-( const Matrix<T>& m ) const
{
    if ( rows != m.rows || cols != m.cols )
        throw SizeMismatch();

    Matrix<T> w( rows, cols );

    for ( int i = 0; i < rows * cols; i++ )
        w.element[i] = element[i] - m.element[i];

    return w;
}

template<class T>
Matrix<T> Matrix<T>::operator-() const
{
    Matrix<T> w( rows, cols );

    for ( int i = 0; i < rows * cols; i++ )
        w.element[i] = -element[i];

    return w;
}

template<class T>
Matrix<T> Matrix<T>::operator*( const Matrix<T>& m ) const
{
    if ( cols != m.rows )
        throw SizeMismatch();

    int ct = 0, cm = 0, cw = 0;
    Matrix<T> w( rows, m.cols );

    for ( int i = 1; i <= rows; i++ ) {
        for ( int j = 1; j <= m.cols; j++ ) {
            T sum =  element[ct] * m.element[cm];
            for (int k = 2; k <= cols; k++) {
                ct++;
                cm += m.cols;
                sum += element[ct] * m.element[cm];
            }
            w.element[cw++] = sum;
            ct -= cols - 1;
            cm = j;
        }
        ct += cols;
        cm = 0;
    }

    return w;
}

template<class T>
Matrix<T>& Matrix<T>::operator+=( const T& x )
{
    for ( int i = 0; i < rows * cols; i++ )
        element[i] += x;

    return *this;
}

template<class T>
void Matrix<T>::Output( std::ostream& out ) const
{
    int k = 0;

    for ( int i = 0; i < rows; i++ ) {
        for ( int j = 0; j < cols; j++ )
            out << element[k++] << "  ";
        out << std::endl;
    }
}

template<class T>
std::ostream& operator<<( std::ostream& out, const Matrix<T>& x )
{
    x.Output( out );

    return out;
}

#endif

